# estimate John Jackson's model using multinomRob
#

options(width=130);

library(multinomRob)

# read John Jackson's replication dataset (PA 2002, Winter)
#
jjdat <- read.table(file="pareplic.asc", header=T);

#> names(jjdat)
# [1] "Voivodship"    "ID"            "UD.KLD"        "SLD"
# [5] "PSL"           "UP"            "Catholic"      "Other"
# [9] "Constant"      "New.Jobs"      "Unemp"         "SOE"
#[13] "Church.Att"    "Schooling"     "Age"           "Farmer"
#[17] "Village"       "Suchoka.UD"    "Pawlak.PSL"    "Kwasn.SLD"
#[21] "Zioqkowska.UP" "Voters"

jjuse <- jjdat;
#total number of voters
jjuse$TotalVotes   <- jjdat$Voters * 1000;
jjuse$SLD      <- jjdat$SLD/100      *jjuse$TotalVotes;
jjuse$PSL      <- jjdat$PSL/100      *jjuse$TotalVotes;
jjuse$UP       <- jjdat$UP/100       *jjuse$TotalVotes;
jjuse$Catholic <- jjdat$Catholic/100 *jjuse$TotalVotes;
jjuse$Other    <- jjdat$Other/100    *jjuse$TotalVotes;
jjuse$UD.KLD   <- jjdat$UD.KLD/100   *jjuse$TotalVotes;

jjuse$Age  <-  jjdat$Age /10;

## Set parameters for Genoud
zz.genoud.parms <- list( pop.size             = 10000,
                        wait.generations      = 30,
                        max.generations       = 100
                        )

# estimate the model using equality constraints, to get multinomT starting values

muljj2 <-
  multinomRob(
    list(
      SLD ~      New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Kwasn.SLD + Suchoka.UD ,
      PSL ~      New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Farmer + Village + Pawlak.PSL + Suchoka.UD ,
      UP ~       New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Zioqkowska.UP + Suchoka.UD ,
      Catholic ~ New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Suchoka.UD ,
      Other ~    New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Suchoka.UD ,
      UD.KLD ~   0 ),
    jjuse,
    equality=list(
    list(
      SLD ~      Suchoka.UD + 0,
      PSL ~      Suchoka.UD + 0 ,
      UP ~       Suchoka.UD + 0 ,
      Catholic ~ Suchoka.UD + 0 ,
      Other ~    Suchoka.UD + 0 )),
    multinom.t = 2,  multinom.t.df = 10000,
#    genoud.parms = zz.genoud.parms,  # use the defaults
    print.level = 3);

summary(muljj2, weights=TRUE);

save(muljj2, file="jj2.RData")
